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Thermal stability properties and the melting-like transition of Najv, N = 13-147, clusters are stud- 
ied through microcanonical molecular dynamics simulations. The metallic bonding in the sodium 
clusters is mimicked by a many-body Gupta potential based on the second moment approximation 
of a tight-binding Hamiltonian. The characteristics of the solid-to-liquid transition in the sodium 
clusters are analyzed by calculating physical quantities like caloric curves, heat capacities, and 
root-mean-square bond length fluctuations using simulation times of several nanoseconds. Distinct 
melting mechanisms are obtained for the sodium clusters in the size range investigated. The calcu- 
lated melting temperatures show an irregular variation with the cluster size, in qualitative agreement 
with recent experimental results. However, the calculated melting point for the Nass cluster is about 
40 % lower than the experimental value. 
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La fusion y las propiedades de estabilidad termica de cumulos de Najv, TV = 13-147, se estudian 
utilizando simulaciones de dinamica molecular en el ensamble microcanonico. El enlace metalico 
en los cumulos de sodio se modela con un potencial de Gupta de muchos cuerpos que se basa en 
la aproximacion de segundo momento de un hamiltoniano de amarre fuerte. Las caracteristicas de 
la transition solido a h'quido en los cumulos de sodio se analizan mediante el calculo de cantidades 
fi'sicas como la curva calorica, el calor espeeffico y la desviacion cuadratica media de las fluctuaciones 
en las distancias interatomicas, utilizando tiempos de simulation de varios nanosegundos. Mecanis- 
mos diferentes de fusion se obtuvieron para cumulos de sodio en el rango de tamanos investigados. 
Las temperaturas de fusion calculadas muestran una variation irregular como funcion del tamano 
del cumulo, en acuerdo cualitativo con resultados experimentales recientes. Sin embargo, el punto 
de fusion calculado para el cumulo de Nass es aproximadamente 40 % mas bajo que el valor exper- 
imental. 

Descriptores: Cumulos metalicos; cumulos de sodio; fusion en cumulos; transition de fase en cumulos 
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I. INTRODUCTION 

The study of the thermal stability and melting transi- 
tion of sodium clusters is nowadays a very active field of 
research. This interest is motivated, in part, by the fact 
that just recently it was possible to make direct compar- 
isons between theoretical results, mainly extracted from 
computer simulations , and experimental data g [l^] 
obtained for the melting-like transition of sodium clus- 
ters of different sizes. In fact, during the last few years 
the caloric curve, heat capacity, and melting tempera- 
tures of Na^v, N = 55-200, clusters have been measured 
using the temperature dependence of the fragmentation 
mass spectra [p~0|— [l2[| . The melting points of these clus- 
ters were found to be on average 33 % (120 K) lower 
than the bulk value, and more surprisingly, large varia- 
tions were observed in the melting temperatures (± 30 
K) with changing cluster size fli"|| . 

On the theoretical side, molecular dynamics (MD) and 
Monte Carlo (MC) methods have been used to provide 
microscopic descriptions on the melting mechanisms of 



sodium clusters in the size range of 8-147 atoms. In these 
studies, the metallic bonding of the sodium clusters has 
been described using different levels of approximation, 
from empirical many-body model potentials and 
tight-binding hamiltonians |3|,|7j to first-principles density 
functional theory @,|j§f|. 

Despite the large amount of information obtained from 
the above theoretical studies, several questions on, for 
example, the irregular variations of the melting temper- 
atures with respect to the cluster size remain unsolved. 
One difficulty existing in the theoretical approaches to 
study the melting-like transition of sodium clusters is re- 
lated with the interplay between the geometric and elec- 
tronic structure effects in such systems |s||],[ll| • Although 
MD simulations that use density functional theory explic- 
itly take into account the electronic degrees of freedom 
and their effects in the cluster geometry, the limitation of 
this approach is related to the relatively short simulation 
times (a few picoseconds) during which time-averages of 
the structural and dynamical properties are calculated. 
This problem, caused by the extremely large computa- 
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tional effort involved in the first-principles MD simula- 
tions, is especially critical in the study of the melting 
transition where large thermal fluctuations are present, 
and therefore much longer simulation times are required 
to calculate converged time averages. 

On the other hand, MD simulations based on many- 
body model potentials allow extension of the simulation 
time up to the nanosecond regime, employing reasonable 
computational resources. However, in this case the de- 
scription of the metallic bonding does not explicitly in- 
clude the electronic degrees of freedom. In the present 
work, we adopt the latter approach to perform molec- 
ular dynamics simulation of the melting-like transition 
of Najv, N = 13, 20, 55, 135, 142, and 147, clusters us- 
ing a many-body Gupta potential 1 13 - 15 1 , and simulation 
times of ~ 50 nanoseconds. The objectives of this study 
are: (1) to test the phenomenological many-body Gupta 
potential in the description of the melting transition us- 
ing adequate simulation times and (2) to compare the 
predictions of this simulation with those obtained from 
the same many-body potential but using the MC method 
for the averaging procedure, and also with the re- 
sults obtained from first-principles MD |6|J^] using shorter 
simulation times. These tests and comparisons will pro- 
vide additional insights on the melting mechanisms in 
sodium clusters and provide useful information on the 
performance of the different simulation methods. 

In Section II we describe the model potential and pro- 
vide details on the simulation method. Results on the 
caloric curves, heat capacities and thermal stability prop- 
erties of different cluster sizes are given in Sec. III. A 
summary of this work is presented in Sec. IV. 



II. METALLIC POTENTIAL AND SIMULATION 
PROCEDURE 

The many-body model potential used in this work cor- 
responds to the Gupta potential that is based on the 
second moment approximation of a tight-binding hamil- 
tonian mM . Its analytical expression is given by: 
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where r , A, £, p, and q are adjustable parameters p4[ . 
For sodium clusters these parameters have been fitted 
to band structure calculations p5[. Their values are: 
A=0.01595 eV, £=0.29113 eV, r =6.99 bohr, p=10!3, 
and 5=1.30 ftq ]. This type of potential has been exten- 
sively utilized in other metal cluster simulations pel poi . 



obtaining good agreement with results generated from 
first-principles methods. 

In order to study the cluster melting-like transition, 
we use the constant-energy MD method to calculate the 
structural changes as a function of the cluster energy. 
Within this method, Newton's equations of motion are 
solved for each atom within the cluster using the Verlet 
algorithm Hj. Through this procedure, one obtains the 
atomic positions and momenta as a function of time, that 
are used to calculate time-averages of physical quantities 
characterizing the cluster structure and dynamics. A typ- 
ical calculation consists in heating up a cluster from its 
zero temperature configuration until it transforms into a 
liquid-like cluster. To simulate this procedure the cluster 
total energy is increased in a step-like manner by scaling 
up the atomic velocities and therefore augmenting the 
kinetic energy of the cluster. 

The simulation starts by slightly perturbing the coor- 
dinates corresponding to the lowest-energy configuration 
of the cluster and setting up the atomic momenta to a 
zero value, in order to eliminate the translation of the 
cluster center of mass and its rotational motion. The 
time step of the MD runs is 2.4 fs, which provides to- 
tal energy conservation within 0.001 %. For each initial 
condition the cluster is equilibrated during 10 time steps 
and the time averages are calculated using 10 7 time steps 
in the solid and liquid regions, but 2 times longer simula- 
tion times are used in the melting-like transition region. 
This amount of averaging time is a necessary condition 
to obtain converged physical quantities characterizing the 
melting process in these finite systems. 

To obtain the lowest-energy structure of each clus- 
ter size we combine simulated quenching techniques |T^ ] 
and evolutive algorithms p2[ , which are able to perform 
global searches of the potential energy surface (PES) in a 
very efficient way despite the large number of degrees of 
freedom involved in the cluster structure optimization. 
These methods do not only provide the lowest-energy 
configuration but the distribution of isomers in the low 
energy range. This optimization procedure has been suc- 
cessfully used in other cluster structure optimizations 

mum- 

The behavior of the structural and thermal properties 
during the heating process of the cluster is monitored 
by calculating the temperature, heat capacity, and root- 
mean-square (rms) bond length fluctuations as a function 
of the cluster total energy using the following expressions: 
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where £/. is the cluster kinetic energy, ks is the Boltz- 
mann constant, <...> denotes a time average, and 
corresponds to the distance between atoms i and j. The 
above mathematical expressions were introduced in Refs. 
p3| , p4j in order to calculate structural and thermal prop- 
erties of clusters from computer simulation studies. 



volume isomerizations have also been obtained for a 13- 
atom cluster modeled by a pair-wise additive potential 
with a soft repulsive core interaction p9|. 
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III. RESULTS AND DISCUSSION 

A. Nais 

The lowest-energy configuration of the Nai3 cluster 
corresponds to an icosahedral structure shown in Fig. 
1(a). This geometry was used to initiate the heating up 
procedure through the MD method. The cluster temper- 
ature (caloric curve) and specific heat as a function of the 
cluster total energy, calculated using Eqs. (3) and (4), are 
displayed in Figs. 2(a) and 3(a), respectively. The change 
in the slope of the caloric curve as well as the existence of 
a maximum in the specific heat are characteristics of the 
solid-to-liquid transition in clusters l23| - [25| . These fea- 
tures are clearly seen in Figs. 2(a) and 3(a), indicating a 
melting-like transition in the Nai3 cluster. This transi- 
tion, which occurs over a broad range of energy involving 
one or more intermediate stages (such as isomerizations, 
coexistence of liquidlike and solidlike phases, partial and 
surface melting etc.), has been widely discussed in ear- 
lier studies of atomic clusters Figure 4(a) shows the 
rms bond-length fluctuation, S, as a function of the clus- 
ter total energy, calculated using Eq. (5). It shows two 
abrupt variations at different energies that are in contrast 
with the analog of the Lindemann criterion p7| for bulk 
melting, where a single abrupt increase in S is observed. 

By performing thermal quenchings from configurations 
generated at different total energies during the MD simu- 
lation, the different melting mechanisms occurring in the 
Nai3 cluster can be investigated. It is found that the 
first abrupt change in S at low energy (temperature) is 
due to cluster isomerization involving only surface atoms, 
whereas the second increase at higher energy (temper- 
ature) corresponds to isomerizations where the central 
(bulk like) atom as well as the surface atoms are involved. 
The onset of the surface and volume isomerizations occur 
at T=149 and 226 K, respectively, whereas the tempera- 
ture corresponding to the maximum of the specific heat 
(which indicates the transition to the liquidlike state) is 
at T=260 K. Similar isomerization mechanisms and tran- 
sition temperature values have been obtained from MC 
simulations of Nai3, using the same many-body model 
potential (71. A similar melting behavior has also been 
obtained for M13 and AI13 through MD simulations using 
the Gupta potential p5 28 1. Nevertheless, this melting 



mechanism is not exclusive of metal clusters interacting 
by a many-body Gupta potential since both surface and 



The double icosahedron with a capped atom over the 
central pentagonal ring displayed in Fig. 1(b) corre- 
sponds to the lowest-energy structure of Na2o- By heat- 
ing up this cluster we obtain the caloric curve and specific 
heat displayed in Figs. 2(b) and 3(b), and the rms bond 
length fluctuation shown in Fig. 4(b). In addition to the 
slight change in the slope of the caloric curve at low en- 
ergies, the specific heat shows a shoulder at low energy 
(temperature), before it reaches its maximum at higher 
energy (temperature). A very abrupt increase in the S 
value is obtained at low energy that corresponds, accord- 
ing to the caloric curve, to a cluster temperature of 57 
K. A second abrupt increase in S at around 157 K is also 
obtained after further heating of Na2o- 

The microscopic features characterizing this melt- 
ing behavior can be extracted by performing periodic 
quenchings using the cluster configurations along the 
MD run at a given cluster energy. The analysis of the 
quenched cluster structures indicate that the first abrupt 
increase in S, that is in correspondence with the shoulder 
in the specific heat at low energies, is related with iso- 
merization transitions where the extra atom incorporates 
into the cluster surface, generating a fluctuating six-atom 
ring in the central region of the cluster structure. This 
isomerization is equivalent to a surface reconstruction of 
the icosahedral surface without diffusion. The second 
jump in S is associated to a surface melting stage where 
the double icosahedron transforms into a more compact 
structure containing two internal atoms. At higher ener- 
gies (temperatures) complete melting is observed, char- 
acterized by diffusive motion where the 6 value levels 
off. The two abrupt increases in S obtained for Na2o are 
equivalent to those described for Nai3, except that in 
the larger cluster the first increase corresponds to an iso- 
merization due to the outer atom lying over the cluster 
surface. These phenomena occurring in Nai3 and Na2o 
before they fully melt, showing atomic diffusion, can be 
denominated as a premelting stage. 

The premelting phenomenon was first obtained in MD 
simulations of M14 and M20 ]30|-|32]| and some Be^, 
A=9,ll,12,14, clusters |53]. Later, other metal clusters 
with magic and non-magic number sizes also reveled a 
premelting stage (see, for example, Refs. ]7|,p5||). The 
two-step melting transition described above was also ob- 
tained using the q-jumping MC method and the Gupta 
potential . A good agreement was found in the pre- 
melting and melting temperatures calculated with this 
approach ^,0. However, higher temperature values, for 
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the two transitions, were found by first-principles orbital- 
free MD on Na2o @, which might be due to the much 
shorter simulation time used in such calculation. 



C. Na 55 



first-principles methods, however intense efforts are cur- 
rently being performed to solve this problem. On the 
other hand, further experimental work is expected in the 
near future that confirms the relatively high value of the 
melting temperature of the Na55 cluster 113] . 



Figure 1(c) shows the lowest-energy configuration of 
Na55 corresponding to the two-shell Mackay icosahedron. 
The caloric curve and specific heat obtained by heating 
up this structure are displayed in Figs. 2(c) and 3(c), 
respectively. Figure 4(c) shows the rms-bond length fluc- 
tuation where an abrupt increase is observed at a total 
energy that, according to the caloric curve, indicates a 
melting temperature of 151 K, whereas the maximum in 
the specific heat is obtained at a slightly higher temper- 
ature of 166 K. In contrast to Nai3 and Na2o, the Nass 
cluster shows a single abrupt change in 5. However, by 
visualizing the dynamical trajectories of each atom in the 
cluster at the temperatures within the transition region, 
we found that the melting process develops in several 
stages. In the first one, the most external layer fluc- 
tuates between an incomplete and complete icosahedral 
surface by expelling and receiving atoms back and forth. 
At higher energies an exchange of atoms between the in- 
termediate and most external layers is obseved, and with 
a further increase in energy, fully melting is found, where 
the central atom also contributes to the cluster diffusive 
motion. This complex melting mechanism has been re- 
cently studied in AI55 and other aluminum clusters by 
introducing dynamical degrees of freedom p4| . 

Similar melting stages have also been obtained by MC 
simulation of Na 55 However, in those calculations 

a slightly higher melting temperature of 175 K was re- 
ported. The first-principles orbital-free MD simulation 
of Nas5 also reported j|] a melting transition at 190 K. 
Again, the smaller melting temperature obtained in the 
present work might be due to the much longer simula- 
tion times we have used in our MD simulations as com- 
pared to those used in Ref. [||. Nevertheless, none of 
the melting temperatures calculated by us and other au- 
thors arc in agreement with the experimental value 
(320 K) reported for Nass H| . This discrepancy between 
theory and experiment has not yet been solved since it 
would require a more detailed modeling of the energy 
landscape of Nass, that includes not only information 
on the basins of attraction of the equilibrium structures 
but also on the topology around saddle points connecting 
the lowest-energy minima. Additionally, in a more de- 
tailed description of the potential energy surface it would 
be possible that the global minimum for Nass does not 
correspond to the icosahedral structure but to an un- 
known special geometry with larger thermal stabiltity 
against melting. At present, it represents a theoretical 
challenge to characterize the potential energy surface of 
systems with such number of degrees of fredom using 
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The global minima of the larger sodium clusters inves- 
tigated in this work correspond to icosahedral structures. 
The three-layer Mackay icosahedron shown in Fig. 1(f), 
is the lowest-energy isomer of Nai47. The lowest-energy 
structures of the Nai35 , Nai42 are incomplete icosahedra 
obtained by removing 12 and 5 vertex atoms, respec- 
tively, from the 147 Mackay icosahedron (see Figs. 1(d) 
and 1(e), respectively). Despite the existence of an in- 
complete surface layer in the lowest-energy structure of 
Nai35, Nai42, they show caloric curves, Figs. 2(d-e), and 
rms bond-length fluctuations, Figs. 4(d-e), very similar 
to those obtained by heating up the complete icosahe- 
dron structure of the 147-atom cluster, Figs. 2(f) and 
4(f). The calculated melting temperatures obtained us- 
ing the Lindemann criterion (6 ~ 0.15) for Nai35, Nai42, 
and Nai47 are 135 K, 190 K, and 171 K, respectively. 
These values are smaller than those obtained from the 
maximum of the specific heat (see Table I). By visualiz- 
ing the atomic coordinates as a function of time at the 
energies where the 6 values change abruptly, it is ob- 
served that the melting is initiated at the cluster surface. 
For the three sizes investigated, the atomic mobility in- 
creases with temperature starting from the most exter- 
nal layer and propagating into the internal layers. This 
stage, known as surface melting [j32f , precedes the com- 
plete cluster melting characterized by the diffusive mo- 
tion of all the atoms in the cluster which is observed at 
temperatures where the specific heat is maximum. 

Our calculated melting temperatures for these larger 
sodium clusters are lower than those obtained by other 
authors using the MC method with the same potential 
and the orbital free MD simulations f§ . This differ- 
ence, as in the smaller cluster sizes, we assume is due to 
the much longer simulation times we have used in our MD 
calculations. On the other hand, although our results 
show that the largest melting temperature corresponds 
to the Nai42 cluster, in agreement with the experimental 
data jLl 12|, the absolute values of our calculated melting 
temperatures are about 30 % lower than the experimen- 
tal values [ fulfil^ . These results indicate that the many- 
body Gupta potential, which does not include electronic 
degrees of freedom, only provide a qualitative description 
of the melting of sodium clusters. 
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IV. SUMMARY 



The melting-like transition of Na^, N = 13, 20, 55, 
135, 142, and 147, has been investigated through mi- 
crocanonical MD simulations using a phenomenological 
many-body Gupta potential. The solid-to-liquid transi- 
tion was studied by calculating caloric curves, rms bond- 
length-fluctuations and specific heats. The indicators of 
the cluster melting correspond to changes in the slope of 
the caloric curve, abrupt increases in the 5 values, and 
the existence of maxima in the specific heats. Table I 
shows the melting temperatures calculated for all clus- 
ter sizes using those criteria. The main features com- 
ing out from these data are: (i ) The melting temper- 
atures calculated from the maxima of the specific heats 
are systematically higher than the values obtained using 
the Lindemann criterion, (ii ) There is an irregular vari- 
ation in the melting temperatures as a function of the 
cluster size, the highest value being the one correspond- 
ing to the Nai42 cluster. These results are in qualitative 
agreement with the experimental data jy],[^. (Hi ) The 
calculated melting temperature for the Nass cluster is 
about 40 % lower than the reported experimental value 
|^2| . (iv ) The melting transition in sodium clusters is 
a complex phenomenon that involves several stages in 
which the system undergoes different structural changes 
(isomerization, premelting and surface melting) before it 
shows a diffusive regime characteristic of the liquidlike 
phase. 

A comparison of the present results with those ob- 
tained using the same many-body potential and the MC 
method , and with those generated from orbital-free 
MD simulations |]||, indicate that in general, the melt- 
ing temperatures calculated by heating up the lowest- 
energy isomer, are lower when much longer simulation 
time is employed. In this work, the simulation time was 
extended up to the nanosecond time regime where it is 
very likely that the time-averages of the physical quanti- 
ties that characterize the cluster melting might be much 
better converged. 

In order to obtain a better (quantitative) agreement 
with the experimental results on the melting of sodium 
clusters it would be necessary to either extend the sim- 
ulation time in the first-principles MD calculations or 
to design a many-body potential that describes with a 
higher level of approximation the complex topology of 
the potential energy landscape of sodium clusters. Work 
in both directions is currently under progress. 
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TABLE I. Binding energies (BE) in eV/atom and melting 
temperatures in K calculated from the temperature value at 
the maximum of the specific heat and using the Lindemann 
criterion (8 ~ 0.15). For N = 13 and 20 there are two values 
due to the existence of two-stage (premelting and melting) 
transitions. 
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0.734 
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0.929 
0.933 
0.935 



260 
220 
166 
181 
189 
180 



149, 226 
57, 157 

151 

135 

190 

171 



320 
250 
285 
272 



a Refs. (lljjlj. 



FIG. 1. Lowest-energy structures of Najv, iV = 13 (a); 20 
(b); 55 (c); 135 (d); 142 (e); and 147 (f); clusters. The clus- 
ter geometries correspond to Mackay icosahedra for N — 13, 
55, and 147; a capped double icosahedron for N = 20; and 
incomplete three-layer Mackay icosahedra for N — 135 and 
142. 



FIG. 2. Caloric curves of Na N , N = 13 (a); 20 (b); 55 (c); 
135 (d); 142 (e); and 147 (f); clusters. The cluster energy 
is calculated taking as reference the value of the binding en- 
ergy of the most-stable (lowest-energy) configuration given in 
Table I. 



FIG. 3. Specific heats of Najv, N = 13 (a); 20 (b); and 
55 (c); clusters. The cluster energy is calculated taking as 
reference the value of the binding energy of the most-stable 
(lowest-energy) configuration given in Table I. 



FIG. 4. RMS bond-length fluctuations of Najv, N = 13 
(a); 20 (b); 55 (c); 135 (d); 142 (e); and 147 (f); clusters. The 
cluster energy is calculated taking as reference the value of 
the binding energy of the most-stable (lowest-energy) config- 
uration given in Table I. 
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